clc; clear; close all;

%% Load data from the microstructure processing routines
load('intermediate\MicrostructureEstimatesAndOutput.mat');

%Run file to input cross-currency basis data
run functions\CIPGraph3MEURCHF;

%Extract from the CIPGraph3MEURCHF processing the various cross-currency
%bases
dateVecCIP = CIPAllS1{:,1};
CIPCalcsAll = [CIPCalcAUD CIPCalcGBP CIPCalcCAD CIPCalcEUR CIPCalcJPY CIPCalcEURCHF];

%Plot the cross-currency bases
plot(dateVecCIP,CIPCalcsAll)
dateaxis('x',2)

CIPVecList = {CIPVecAUD,CIPVecGBP,CIPVecCAD,CIPVecEUR,CIPVecJPY,CIPVecEURCHF};
counter = 1;
dateVecCIP = dateVecCIP(3077:end,:);
dateVecCIP = dateVecCIP(selectorNonHolidays,:);

%Initialize matrix to store welfare and arb size calculations, matrix
%alternates between welfare and arb size computations, moving through each
%market (column 1 = AUD/USD CIP Welfare, column 2 = AUD/USD CIP Arb
%size...)
WelfareArbs = zeros(305,12);
FunctionEvals = zeros(305,6);

%Initialize vectors with contract parameters given by contract
%specifications, including tick size (minimum fluctuation) and size of
%contracts (how many dollars a contract is worth typically during our estimation period)
TickSizes = [0.00005 0.0001 0.00005 0.00005 0.0000005 0.0001];
ContractSizes = [75000 81250 75000 150000 120000 150000];

time = 0;
%Cycle through each contract and compute the welfare gains to arbitrage and
%gap-closing size
while(counter<=6)
    %Extract the interest rates, spot, and forward data
    CIPVec = CIPVecList{1,counter};
    %Isolate only the recent period for which we have microstructure
    %estimates by day
    CIPVec = CIPVec(3077:end,:);
    %Create a selector vector that eliminates holidays and select out data
    selectorCIP = zeros(315,1)+1;
    selectorCIP([8 9 13 65 85 130 195 260 270 275],:) = 0;
    selectorCIP = (selectorCIP>0);
    CIPVec = CIPVec(selectorCIP,:);
    
    time = 0;
    %Cycle through each day in the recent microstructure dataset and
    %compute welfare gains to arbitrage and size of trade
    while(time<size(CIPVec,1))
        time = time + 1;
        %One can try varying the parameter to calculate welfare since this
        %is the cross-price impact assumption: robust=10 implies the price
        %impact is only 1/10 of what it normally would be due to the fact
        %that spot and forward markets feature significant spillove price
        %impact (going short one leg drives the other leg down and vice
        %versa, resulting in largely cancelling effects from the trade)
        robust = 10;
        robust2 = 100;
        %Calculate assuming (zero coefficient) no interest rate price
        %impact for simplicity since it is second order effect
        Values = @(m) [(CIPVec(time,1)+sign(m)*(TickSizes(1,counter)*coeffs(time,counter)/robust).*(abs(m)*1000000/(ContractSizes(1,counter))).^0.5)./(CIPVec(time,2)-sign(m)*(TickSizes(1,counter)*coeffs(time,counter)/robust).*(abs(m)*1000000/ContractSizes(1,counter)).^0.5) .* (1./(1+(CIPVec(time,4)/100 - sign(m)*(0.0000000/robust2).*(abs(m)*1000000/120000).^0.5733))).^0.25 - (1./(1+(CIPVec(time,3)/100+sign(m)*(0.0000000/robust2).*(abs(m)*1000000/100000).^0.5852))).^0.25]';
        %Values2 is not used
        Values2 = @(m) [(CIPVec(time,1)+sign(m)*(TickSizes(1,counter)*coeffs(time,counter)/robust).*(abs(m)*1000000/(ContractSizes(1,counter))).^0.5) (CIPVec(time,2)-sign(m)*(TickSizes(1,counter)*coeffs(time,counter)/robust).*(abs(m)*1000000/ContractSizes(1,counter)).^0.5)  (CIPVec(time,4)/100 - sign(m)*(0.000000013/robust2).*(abs(m)*1000000/120000).^0.5733) (CIPVec(time,3)/100+sign(m)*(0.000000/robust2).*(abs(m)*1000000/100000).^0.5852)]';
        %Set optimization options
        options = optimoptions('fsolve','FunctionTolerance',1e-12,'OptimalityTolerance',1e-24,'StepTolerance',1e-12,'MaxFunctionEvaluations',1e3);
        %Use solver to find gap-closing size
        [m1,m1Fn] = fsolve(Values,1,options);
        %Numerical integration routine
        grid = (1:1:abs(m1))*sign(m1);
        evals = Values(grid);
        TotalWelfare = sum(evals);
        %Store computations
        WelfareArbs(time,1+(counter-1)*2)= TotalWelfare;
        WelfareArbs(time,2+(counter-1)*2)= m1;
        FunctionEvals(time,counter)= m1Fn;
    end
    counter = counter + 1;
end

mean3Ms = mean(CIPCalcsAll);
std3Ms = std(CIPCalcsAll);
min3Ms = min(CIPCalcsAll);
max3Ms = max(CIPCalcsAll);
meanabs3Ms = mean(abs(CIPCalcsAll));

%These are the data used for the CIP Calculations:
%ADBB3M CMPN Curncy Australian Bank Bill 3 Month
%BP0003M Index STERLING ICE LIBOR
%CDOR03 Index 3 Month Canadian Banker’s Acceptance Rate
%GTDEM3M Corp German Treasury Bill (Bubill)
%GB03 Govt (Three month T-Bills)
%JY0003M Index Japan 3 Month LIBOR


data2 = [mean3Ms;std3Ms;min3Ms;max3Ms;meanabs3Ms];
input.data = data2;
input.tablePositioning = 'h';
input.tableColLabels = {'AUD/USD','GBP/USD','CAD/USD','EUR/USD','JPY/USD','EUR/CHF'};
input.tableRowLabels = {'Mean 3M Cross Currency Basis','Std 3M Cross Currency Basis','Min 3M Cross Currency Basis','Max 3M Cross Currency Basis','Mean Abs 3M Cross Currency Basis'};
input.dataFormatMode = 'column';
input.dataFormat = {'%.4f',5};

counter = 1;
WelfareArbs2 = zeros(3391,12);
FunctionEvals2 = zeros(3391,6);
%Cycle through each contract and compute the welfare gains to arbitrage and
%gap-closing size
time = 0;
while(counter<=6)
    CIPVec = CIPVecList{1,counter};
    %Non-crisis from Du paper is t=481 (01/01/2010) on to t=3088 (12/31/2019)
    time = 480;
    %This utilizes the full dataset as opposed to recent data only in
    %earlier loop
    while(time<3088)
        time = time + 1;
        %One can try varying the parameter to calculate welfare since this
        %is the cross-price impact assumption: robust=10 implies the price
        %impact is only 1/10 of what it normally would be due to the fact
        %that spot and forward markets feature significant spillove price
        %impact (going short one leg drives the other leg down and vice
        %versa, resulting in largely cancelling effects from the trade)
        robust = 10;
        robust2 = 10;
        %Calculate assuming (zero coefficient) no interest rate price
        %impact for simplicity since it is second order effect
        Values = @(m) [(CIPVec(time,1)+sign(m)*(TickSizes(1,counter)*alphas(1,counter)/robust).*(abs(m)*1000000/(ContractSizes(1,counter))).^0.5)./(CIPVec(time,2)-sign(m)*(TickSizes(1,counter)*alphas(1,counter)/robust).*(abs(m)*1000000/ContractSizes(1,counter)).^0.5) .* (1./(1+(CIPVec(time,4)/100 - sign(m)*(0.0000000/robust2).*(abs(m)*1000000/120000).^0.5733))).^0.25 - (1./(1+(CIPVec(time,3)/100+sign(m)*(0.0000000/robust2).*(abs(m)*1000000/100000).^0.5852))).^0.25]';
        %Values2 not used
        Values2 = @(m) [(CIPVec(time,1)+sign(m)*(TickSizes(1,counter)*alphas(1,counter)/robust).*(abs(m)*1000000/(ContractSizes(1,counter))).^0.5) (CIPVec(time,2)-sign(m)*(TickSizes(1,counter)*alphas(1,counter)/robust).*(abs(m)*1000000/ContractSizes(1,counter)).^0.5)  (CIPVec(time,4)/100 - sign(m)*(0.000000013/robust2).*(abs(m)*1000000/120000).^0.5733) (CIPVec(time,3)/100+sign(m)*(0.00000021/robust2).*(abs(m)*1000000/100000).^0.5852)]';
        %Set optimization options
        options = optimoptions('fsolve','FunctionTolerance',1e-12,'OptimalityTolerance',1e-24,'StepTolerance',1e-12,'MaxFunctionEvaluations',1e3);
        %Use solver to find gap-closing size
        [m1,m1Fn] = fsolve(Values,1,options);
        %Numerical integration routine
        grid = (1:1:abs(m1))*sign(m1);
        evals = Values(grid);
        TotalWelfare = sum(evals);
        %Store calculations
        WelfareArbs2(time,1+(counter-1)*2)= TotalWelfare;
        WelfareArbs2(time,2+(counter-1)*2)= m1;
        FunctionEvals2(time,counter)= m1Fn;
    end
    counter = counter + 1;
end


%Plot values
dateVecCIP2 = CIPAllS1{:,1};
plot(CIPAllS1{481:3088,1},WelfareArbs2(481:3088,[1 3 5 7 9 11]));
dateaxis('x',2)

%New Table Computation of Summary Statistics for Microstructure Dataset
TickSizes = [0.00005 0.0001 0.00005 0.00005 0.0000005 0.0001];
Multipliers = [100000 62500 100000 125000 12500000 125000];
DivideFactor = [20000 10000 20000 20000 2000000 20000];
TickValue = [];
valsForSumStats2 = zeros(8,5);
counter = 0;
while(counter<5)
    counter = counter+1; 
    AdditionalValues = AdditionalValuesAll{1,counter};
    innercounter = 0;
    daVector = zeros(315,6);
    daPointer = 1;
    while(innercounter<5)
        innercounter=innercounter+1;
        daVector(daPointer:daPointer+size(AdditionalValues{innercounter,1},1)-1,:)=AdditionalValues{innercounter,1};
        daPointer=daPointer+size(AdditionalValues{innercounter,1},1);
    end
    daVector=daVector(selectorNonHolidays,:);
    %Format is minimum, mean, maximum price, then daily $ volume, then min,
    %mean, and max volume, then mean number of transactions
    valsForSumStats2(:,counter)=[min(daVector(:,1))/DivideFactor(counter) mean(daVector(:,2))/DivideFactor(counter) max(daVector(:,3))/DivideFactor(counter) mean(daVector(:,4))*Multipliers(counter)*TickSizes(counter) prctile(daVector(:,5),10) prctile(daVector(:,5),50) prctile(daVector(:,5),90) mean(daVector(:,6))]';
end

%Plot curves that illustrate price impacts of different trade sizes in the
%various currency markets
PriceImpactFunctions = zeros(1000,5);
Levels = (1e8:1e8:1e11);
counter = 0;
%Cycle through all markets currently in paper, which is everything except
%CHF/USD
while(counter<5)
   counter = counter+1;
   contracts = Levels./(Multipliers(counter)*valsForSumStats2(2,counter));
   PriceImpactFunctions(:,counter) = (100*TickSizes(counter)/valsForSumStats2(2,counter))*(alphas(1,counter)*(contracts.^0.5))';
end

save('intermediate\CompileAllCIPCalcs');
